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We show that the equations underlying the GW approximation have a large number of solutions. 
This raises the question: which is the physical solution? We provide two theorems which explain 
why the methods currently in use do, in fact, find the correct solution. These theorems are general 
enough to cover a large class of similar algorithms. An efficient algorithm for including self-consistent 
vertex corrections well beyond GW is also described and further used in numerical validation of the 
two theorems. 



INTRODUCTION 

The GW approximation [1 is a many-body technique 
used typically for the calculation of spectral density func- 
tions of solids. Its accuracy for the determination of band 
gaps of insulators as well as its parameter-free nature 
makes it a very attractive method in condensed matter 
physics. The self-consistent GW approximation is a fixed 
point method involving multidimensional objects like the 
Green's function and the self-energy, and it was recently 
demonstrated for an artificial one-point-model that this 
fixed point is not unique and that a different way of iter- 
ating the equations leads to a different solution [2]. It was 
further argued that including vertex corrections could ex- 
acerbate this non-uniqueness problem and lead to several 
solutions. Given that the GW approximation is a state 
of the art method for band structure calculations, such 
an ambiguity is a serious issue. It is also important to 
understand how one can go beyond GW without running 
into unphysical solutions. 

In this article we describe a new algorithm for comput- 
ing the self-energy well beyond the GW approximation. 
We investigate the nature of the additional solutions nu- 
merically and further provide general theorems, that ex- 
plain why the methods currently in use to solve the GW 
equations do indeed lead to a unique solution. Since we 
cover a large class of approximations these results not 
only validate the GW calculations that are done, but they 
also provide conditions on approximations going beyond 
GW for obtaining a meaningful result. 

The starting point are the Hedin equations, which ap- 
pear as Eqs. (A22)-(A25) in the appendix of the 1965 
article of Hedin[3]. We rewrite them here in modern no- 
tation: 
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In these equations G is the Green's function, W is 
the renormalized Coulomb propagator @J, E is the self- 
energy, II is the polarization, r (l,2;3) = A5(l, 2)5(1, 3) 
is the bare vertex, T is the renormalized vertex, A is the 
coupling constant, and the potential differential 5V is the 
sum of the external and Hartree contributions. The no- 
tation (1) =x% = (ri,a\,tx) is used throughout [5]- Note 
that we introduced the coupling constant A such that in 
the non-interacting case the vertex function vanishes [5]. 
Instead one could also define it in such a way that the 
Coulomb propagator vanishes. Using simple transforma- 
tions (introduced later in Eq. (12); use a = 1, b — 
A 2 , c = A -1 ) one can show that these definitions are 
in fact equivalent. The physical meaning of this is, that 
it does not matter whether we define the non-interacting 
limit by switching off the interaction of the electrons with 
the photons or by setting the photon propagator to zero. 

Hedin has shown, how one can gain an expansion of the 
vertex and hence of E and II in terms of the renormalized 
quantities G and W using these equations. That way one 
gets E and II as functionals E[G, W] and In 
addition to the five Hedin equations are the two coupled 
Dyson equations: 

G(1,2) = G (1,2)+ / G (l,3)E(3,4)G(4,2)d(3)d(4) 



W(1,2) = W (1,2) 



Wb(l, 3)11(3, 4)W(4, 2)d(3)d(4), 
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where Go is the Green's function of the non-interacting 
system (which includes the Hartree potential) and 
Wo(l,2) = t^£~t is the bare Coulomb propagator. 
Solving the Dyson equations in conjunction with the 
Hedin equations yields the functionals G[Gq,Wo] aid 
W[G ,W }. 

The separation of the problem into equations for E and 
n as functionals of G and W, as well as G and W as func- 
tionals of Go and Wq is an important conceptual step. In 
a later article by Hedin and Lundqvist [6] , the equations 
are combined and the functional derivative 5E/5G is in- 
troduced. We would like to stress that this vertex equa- 
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tion together with Eqs. |2]), ([3| and the Dyson equations 
are not immediately useful. One also needs the equations 
Q and ([5]) in order to get an expansion of the vertex be- 
yond GW. 

ALGORITHMS FOR HEDIN'S EQUATIONS 

Almost all practical calculations of Hcdin's equations 
use the GW approximation. This amounts to approx- 
imating the full vertex T by the bare vertex, and thus 
the self-energy takes on the simple form £(1,2) — 
i\ 2 G(l, 2)W(1, 2). In the present work we describe a new 
algorithm for solving Hedin's equations which includes 
corrections far beyond the GW approximation. From 
the outset we will insist that the computational storage 
requirements scale as N 3 and the number of operations 
scale as N 4 . These conditions are met by Algorithm [l] 
(see structogram) for calculating E[G, W] and H[G, W] 
using a non-trivial vertex. Note that all the relations in 
the algorithm are exact apart from the two which have 
the derivatives of T removed - these would require TV 
storage. Solving Algorithm [I] together with the Dyson 
equations we refer to as the 'Starfish' algorithm. It is 
straight-forward to work out which diagrams this algo- 
rithm corresponds to: finding the the self-consistent so- 
lution to Starfish is equivalent to solving 




for the vertex, using this to find S and II and subse- 
quently G and W. The entire procedure is performed 
sclf-consistcntly. 

NUMBER AND STABILITY OF SOLUTIONS 

Discretization of space-time is required for the purpose 
of closely examining solutions to these equations. There 
is some ambiguity as to how this should be done, but at 
this level, we merely assume that there are N space-time 
points in total, and that limits in the time variable such 



Algorithm 1 Hedin equations solver for £[G, W\ and 
n[G,W]. Here the shorthand 
G'(l,2;3) = 6G(12)/5V(3), etc., is used and the 
function arguments arc omitted. 
Require: G and W 

Set r = r 

repeat 

G' = J GGV 

W = -iXf (CGF + GG'T + fi/fiffff) 

w' = / www 

£' = iXf (w'GT + WGT+ W/ffl/Tf) 
r = r + £' 

until r converged 
£ = iX JWGT 
n = —iX J GGF 



as G(4, 2 + ) are taken to mean G(4,2). Also the Dirac 
delta function in Eq. ([I]) now becomes a Kronenker delta 
function. In doing so one loses the physical meaning of 
these equations, but we will assume that the true physical 
solution can be recovered in the continuum limit when 
N -> 00. 

Irrespective of whether the GW approximation, 
Starfish or some other truncation is used, the equations 
to be solved form a closed system of polynomial equa- 
tions enabling us to prove general theorems about such 
a system. 

To do so we define a concise notation. Let 

f= (G,w,n,E,r) e c n , 

be the vector of all dependent (or unknown) variables. 
Here n = 47V 2 + N 3 is the number of unknowns. We 
regard the tensors Go and Wq (and the bare vertex) as 
known and fixed. Equations pi), ([3]), ^ and a vertex 
equation like Eq. ^ can now be written in a compact 
form as F = q(F) or t)(F) = with t)(F) = g(F) - F, 
where q = (qi, . . . , g n ) is a set of polynomials in several 
variables. 

It is important to point out that most of the following 
considerations do not depend on the precise form of the 
vertex equation. For example the trivial vertex equation 
T = T corresponding to the GW approximation is also 
allowed here. In this case one can also eliminate the ver- 
tex from the equations and redefine F, q and t) in order to 
include only the smaller set of quantities and equations. 
Similarly, if one were to fix W in the Starfish algorithm 
then F would be (G, E, T) and the equations for W and 
n could be eliminated. 

Since we are interested in the dependence of these 
equations and their solutions on the coupling strength 
A, the equations to be solved are 

F = q x {F) or t) X (F) = 0. (8) 
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Multiple fixed points 



We now want to determine the number of solutions 
to Eq. (pi). An upper bound is provided by Bezout's 
theorem]?] which states that the maximum number of 
solutions to a system of polynomial equations (if finite) 
is equal to the product of the total degree of each equa- 
tion. Thus the GW approximation with fixed W has at 
most 2 N solutions, and the Starfish algorithm, also for 
fixed W, has at most 7^ 2 2N solutions. The Bezout 
bound fails to take into account sparsity in the poly- 
nomial equations and therefore also degeneracy of the 
roots. Buchberger's algorithm is a method of system- 
atically determining the exact number of roots by de- 
composing the equations into a Grobner basis [8]. This 
procedure is, however, computationally very demanding 
and can be performed only for small (and therefore non- 
phyiscal) N. For example, when N — 2 the GW approx- 
imation with fixed W has precisely 6 solutions. Likewise, 
for N — 1 the Starfish algorithm yields 3 solutions. From 
these considerations, it seems quite surprising that self- 
consistent GW works at all for realistic values of N. We 
will now provide two theorems that may explain this ap- 
parent success. 



Theorem 1. Let f) ; 

polynomials fi£ : C" 



= (f, A ,...,f) A ) 7 AG [0,1] be a set of 
C with the following properties: 



(i) They depend pointwise continuously on X. 

(ii) For vanishing A there is exactly one solution F° to 
t)°(F)=0. 

(Hi) The Jacobian J A -(F) = §^ satisfies det[J°(F )] ^ 
0. 

Then we have 

(a) For every ball B r centered at Fq there is a Ao > 
such that for all smaller X, t) A also has a zero in 
B r . 

(b) For every ball Br centered at Fq there is a Ao > 
such that for all smaller X, f) A has at most one zero 
in B R . 

It can be checked easily, that these conditions are sat- 
isfied by all versions of the f) A we defined above. The 
situation is depicted in Fig. [I] If the interaction is small, 
there is only one solution -Fphys that is close to the non- 
interacting one. All others tend to infinity in the limit of 
small interaction, i.e. for A — > one can let the radius r 
tend to zero and R tend to infinity. This behavior sug- 
gests, that at least in some low coupling regime -Fphys is 
indeed the physical solution, while all other fixed points 
are far away from the right result. 

Proof, (a) First we note, that owing to the continuity 
with respect to A given in (i) and the fact that t) A is 




FIG. 1: Sketch of the position of the solutions of the 
interacting system relative to the non-interacting one 
F°. For small interaction there is one solution 
-Pphys € B r of Eq. ([8| , that is close to the 
non-interacting one F°. All others are outside Br, 
while r tends to zero and R tends to infinity for the 
small coupling limit A — > 0. 



analytic then J X (F) also depends continuously on both 
A and F. So using (iii), we can restrict F and A to 
be in a region around F and 0, respectively, such that 
J A (F) is close enough to J°(F°) to have non vanishing 
determinant 



det[J A (i^)] ^ 0. 

We can assume, that the ball B r named in (a) was already 
small enough to ensure this. Now we define 



min {||t)°(F)||} ^0, 

FedB r 1 ' 



where dB r denotes the surface of B r and ||-|| denotes 
the Euclidean vector norm. Now we are interested in 
the minimum of ||f) A || on B r = B r U dB r . Clearly this 
minimum tends to zero for A — > since t) a (F ) — 0. On 
the other hand, the minimum of ||f) A || with respect to F 
on the surface dB r tends to m > 0. So we can further 
restrict A such that ||f) A || has at least one local minimum 
in the interior of B r . Let's call the location of one such 
minimum -Fphyg. Now that we know that this minimum 
is not on the surface of B r , but in the interior, we can 
conclude 



d\\t) 



A 1 1 2 



d(Re[Fi\) 



0. 



d\\t) 
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d(lm[Fi]) 



0. 
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This implies 
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X! Jji(Fphys)[t)j(F p hys) 



= 0. 



(9) 
(10) 



Now we have chosen A and the size of the ball B r such 
that the Jacobian has non-zero determinant. So the only 



way Eq. ( 10 ) can hold is if 

f)Vphys) = 0. 



□ 



Proof, (b) We will do this proof in two steps. First we 
will chose an r > such that we can show that a ball with 
this radius, centered at F , contains no more than one 
solution (given A is small enough). Then we will show, 
that for the same ball B r and sufficiently small A there 
is no solution in the set Br \ B r , where Br is the ball 
named in (b). 
Step 1. Define 



M x [q. 



J X {q r (s))ds, 



with q r : [0,1] -> B r . Now M x [q r ] can be decomposed 
into J°(F°) plus a reminder, that vanishes in the limit of 
A — > and r —> 0. That allows us to fix r and restrict A 
such that det(M x [q r ]) ^ for any map q r . Now assume 
we have two zeros F , F 2 £ B r of [) A . Then 



= f) A (F 2 
- 1 d 



f 1 ^-t) X [ S F 2 + (l-s)F 1 }ds 
Jo ds 

£(F 2 - F^jMfrqr], 



where q r (s) = sF 2 + (l — s)F 1 . Now we use that det [M] ^ 
and therefore F 2 = F 1 . So there is indeed no more than 
one solution in B r . 

Step 2. Due to the continuity with respect to A granted 
in (i) the infimum of ||f) A | with F restricted to Br \ B r 
tends to the infimum of ||rj°| with the same restriction 
on F. This is not zero, since we assumed that i)°(F) = 
has only one solution F°. □ 



A natural starting point for this is the non-interacting so- 
lution F — F °. The hope is that this sequence converges 
to a fixed point, i.e. a solution to the equations. A priori 
it is unknown if the calculation will actually converge, or 
if a fixed point obtained this way actually corresponds to 
the physical solution. So one may wonder, why GW and 
similar schemes do, in fact, work in many situations. We 
now show that, for weak coupling once again, a unique 
and convergent solution is guaranteed. 

Theorem 2. If X is not too large, iterating the equations 
according to (11), starting from the initial guess Fq = F° , 
converges to the physical solution -Fphys- 

The connection to the physical solution is made 
through Theorem [T] 

Proof. For a,b,c £M. we define the following transforma- 
tion: 



G' Q = aG W[, = bW Q T' Q = cT 

G' = aG W' = bW r' = c r 
£' = a _1 E n' = 6 _1 n A' = cA 



(12) 



with a 2 be 2 = 1. (We will use these transformations in 
this section only, to avoid confusion of the meaning of 
the primes with the derivatives as used earlier). This 
transformation leaves Eqs. |2|, Q3J) , ^ and e.g. §7§ 
invariant: 

Fl +1 = Q '(F>). 

We now apply the transformation with, say, a — A 1 / 4 , b — 
A 1 / 2 , c = A -1 / 2 . This way g' depends on A explicitly and 
implicitly through a, b and c. Observe, that all coeffi- 
cients appearing in g' tend to zero as A — > 0. The same is 
true for the transformed starting values Fq — F° . In this 
situation Banach's fixed point theorem can be applied for 
the map g' defined in a vicinity of the non-interacting so- 
lution F° . We conclude, that for small A the transformed 
quantities tend to a fixed point. Since the transforma- 
tion can be inverted, this remains true for the original 
quantities. By Theorem [T] for small A the solution -Fphys 
is the solution nearest to the non-interacting one. Hence 
the fixed point obtained is indeed F p i iys . □ 



Generalizations of the Theorems 



Convergence of iterative solutions 

In practice, one solves Eq. Q iteratively. That is, one 
starts with some initial guess Fq, inserts it into the right 
hand side of F = q(F) and obtains a new guess. Iterating 
this procedure defines a sequence 

F i+1 =Q(Fi). (11) 



For the sake of clarity, we presented the Theorems [T] 
and [2] slightly less general than possible. However, these 
Theorems can be easily generalized. For example in the 
proof of Theorem [l] there was no reference to the fact 
that we are dealing with polynomials, rather only to the 
preposition that the functions are analytic in the region 
of interest (i.e. B r and Br). A consequence of this is 
that Theorem [l] also applies to vertex equations that are 
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more sophisticated than just polynomials. It also allows 
us to use Dyson's equations in the 'solved' form, e.g. 
G = (1— GoS) _1 Go, if we restrict E and II (and hence Br 
and B r ) such that the introduced poles are outside Br 
and B r . This restriction is natural in the non-interacting 
limit since then the self-energy and polarization should 
tend to zero. Similarly, one can extend the proof of The- 
orem [2] to these solved Dyson equations (note that they 
still obey the symmetry used in the proof). Also the re- 
striction of choosing the non-interacting solution as start- 
ing point is actually not necessary, though it renders the 
proof much simpler; it can be shown that the fixed point 
becomes attractive in the small coupling limit, even with- 
out the restriction of starting from zero self-energy and 
polarization. 

NUMERICAL INVESTIGATION 

Numerical checks of the above theorems were per- 
formed for both self-consistent GW and the Starfish al- 
gorithm. As already mentioned, obtaining all solutions is 
in practice only possible for very small N . Hence for GW 
we set N — 2 and for Starfish we use N = 1. Plotted 
in Fig. [2] are all the solutions for these algorithms, as a 
function of A. The numerical input, in this case Go and 
W , were chosen to be random complex numbers, and W 
was kept fixed (which is common practice for real GW 
calculations). As mentioned earlier, there are 6 solutions 
in the GW case. Of these, 5 tend to infinity and the 
remaining solution tends to Go as A — > 0. This is a visu- 
alization of Theorem [T] For Starfish, 2 of the 3 solutions 
tend to a constant. This may seem to be in violation of 
the theorem, but in this case the vertex T (and therefore 
F) diverges. 

We can also examine the domains of convergence for 
both of these algorithms. For N — 1 we fix Wo = 1 
and A = 1 and plot the region of convergence of G for 
the fully self-consistent GW and Starfish algorithm (for 
which W is also computed self-consistently) in Fig. [3] It 
can be observed that the region of stability shrinks for 
the higher-order method. Also noteworthy is that the 
region has a fractal boundary (this may be unsurprising 
since for case GqW = 1 the domain is the Mandelbrot 
set). Perhaps more interesting is the region of starting 
points for which the algorithms converge. These are plot- 
ted in Fig. [4] for the same Wo and A but this time with 
G = 1 + i for GW and G = 1/4 + i/4 for Starfish, 
and with a variable starting point for G. Once again 
the region of convergence is smaller for Starfish, but in 
both cases only one solution is found, irrespective of the 
starting point. This is a numerical confirmation of The- 
orem [2] Note that for GW a situation was picked, where 
the non-interacting starting point does not lead to con- 
vergence. Hence this can be considered a large coupling 
situation. But still there seems to be only one stable 




FIG. 2: Plot of the distance of a matrix element of the 
Green's function to the non-interacting one versus the 
coupling strength A for all possible solutions of GW 
with N = 2 and Starfish with N = 1, for random Go 
and W. W is kept fixed in both cases. Always only one 
solution tends to the non-interacting one for the weak 
coupling limit. 
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FIG. 3: Domain of convergence of GW and Starfish 
with N = 1 for input values of Go in the complex plane, 
when using the non-interacting solution as starting 
point. Here Wo = 1 and A = 1. The crosses mark the 
chosen Go for investigating the starting point 
dependence while fixing G , see Fig. |4| 



fixed point. The boundary of the region is also fractal 
(this corresponds to the Julia set). 
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FIG. 4: Domain of convergence of GW and Starfish 
with N — 1 for different starting points of the fixed 
point cycle. The values of Go are fixed to 1 + i for GW 
and 1/4 + i/4 for Starfish, as indicated by the crosses in 
Fig. [3| 

CONCLUSIONS 

We have argued that truncating Hedin's equations to 
some order yields systems of polynomial equations which 
have a very large number of solutions. As an example 
of this, the Starfish algorithm was introduced which in- 
cludes vertex corrections beyond GW and consequently 
has even more fixed point solutions. The number of so- 
lutions tends to infinity as either the order of truncation 
or N tends to infinity, reflecting the inherent problem 
of solving Hedin's equations as a functional differential 
equation. Two theorems were presented that shed some 



light on the general behavior of these fixed points. In par- 
ticular we have shown, that there is exactly one solution 
that tends to the non-interacting case for small coupling, 
while all others are divergent in this limit. Numerical 
tests of self-consistent GW and the Starfish algorithm 
for small N demonstrated that the system also converges 
uniquely to one fixed point even for fairly large coupling. 
Furthermore, the region of stability may be fractal in 
nature, indicating that finding simple necessary and suf- 
ficient conditions for ensuring convergence of GW calcu- 
lations a priori, may be impossible. 

We thank Lucia Reining, Martin Stankovski and Ralph 
Tandetzky for useful discussions. 
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